
##########################################
### OECD regional inequality: analysis ###
##########################################

rm(list = ls())
gc()

library(ggplot2)
library(data.table)
library(boot)
library(reshape2)
library(AER)
library(stargazer)


### Set WD
setwd("~/Dropbox (Princeton)/Research/*Spatial_Inequality/paper_SMD_inequality/JoP_final/replication")  


### -------
### Load OECD data
### -------

load("data/oecd_rwb.RData")


### -------
### CPDS data
### -------

load("data/cpds.RData")

cpds[ , pres_bin := ifelse(pres %in% c("Semi-presidential dominated by president","Presidential system"),1,0)]
cpds[ gov_type %in% c("Single-party majority government","Single-party minority government"), coalition_gov := 0 ]
cpds[ gov_type %in% c("Minimal winning coalition","Surplus coalition","Multi-party minority government"), coalition_gov := 1 ]

cpds <- cpds[ ,.(ISO3,year,fed,pres_bin,bic,eu,gov_left1,coalition_gov)]

### Merge
oecd_rwb <- merge(oecd_rwb, cpds, by=c("ISO3","year"), all.x=T)  

### Full set
oecd_rwb <- na.omit(oecd_rwb)


### -------
### Summary stats 
### -------

d_sum_stat <- unique(oecd_rwb[ ,.(ISO3, country, REG_ID, year, elect_rules, sd_gini_before_tax, diff_gini_idd, year_elec, level_name) ])

sum_stat1 <- d_sum_stat[ ,.(country,elect_rules,sd_gini_before_tax,diff_gini_idd,year_elec,level_name)]
sum_stat2 <- d_sum_stat[ , .N,  country ]

sum_stat_comb <- merge(sum_stat2, sum_stat1, by="country")

sum_stat_comb[ , level_names2 := paste(level_name, " (", N, ")", sep="")]
sum_stat_comb[ , sd_gini_before_tax := round(sd_gini_before_tax, 3)]
sum_stat_comb[ , diff_gini_idd := round(diff_gini_idd, 3)]

sum_stat_comb_fin <- unique(sum_stat_comb[ ,.(country,level_names2,elect_rules,sd_gini_before_tax,diff_gini_idd,year_elec)])


### SI Table A.1
print(xtable::xtable(sum_stat_comb_fin, include.colnames=T, caption="Summary statistics OECD"), include.rownames=F)


sum_stat_covars <- as.data.table(rbind(cbind("var" = "GDP (log)", "mean" = mean(oecd_rwb$gdp_log), "sd" = sd(oecd_rwb$gdp_log), "min" = min(oecd_rwb$gdp_log), "max" = max(oecd_rwb$gdp_log), "source" = "OECD National Accounts"),
                                       cbind("var" = "Unemployment rate", "mean" = mean(oecd_rwb$unemp_rate), "sd" = sd(oecd_rwb$unemp_rate), "min" = min(oecd_rwb$unemp_rate), "max" = max(oecd_rwb$unemp_rate), "source" = "OECD National Accounts"),
                                       cbind("var" = "Trade union density", "mean" = mean(oecd_rwb$trade_union_density), "sd" = sd(oecd_rwb$trade_union_density), "min" = min(oecd_rwb$trade_union_density), 
                                             "max" = max(oecd_rwb$trade_union_density), "source" = "OECD National Accounts"),
                                       cbind("var" = "Pre-tax gini", "mean" = mean(oecd_rwb$gini_market_idd), "sd" = sd(oecd_rwb$gini_market_idd), "min" = min(oecd_rwb$gini_market_idd), "max" = max(oecd_rwb$gini_market_idd), "source" = "OECD IDD"),
                                       cbind("var" = "Post-tax gini", "mean" = mean(oecd_rwb$gini_disp_idd), "sd" = sd(oecd_rwb$gini_disp_idd), "min" = min(oecd_rwb$gini_disp_idd), "max" = max(oecd_rwb$gini_disp_idd), "source" = "OECD IDD"), 
                                       cbind("var" = "Seat share left governments", "mean" = mean(oecd_rwb$gov_left1), "sd" = sd(oecd_rwb$gov_left1), "min" = min(oecd_rwb$gov_left1), "max" = max(oecd_rwb$gov_left1), "source" = "CPDS"),
                                       cbind("var" = "Presidential system (0/1)", "mean" = mean(oecd_rwb$pres_bin), "sd" = sd(oecd_rwb$pres_bin), "min" = min(oecd_rwb$pres_bin), "max" = max(oecd_rwb$pres_bin), "source" = "CPDS"),
                                       cbind("var" = "Bicameralism index", "mean" = mean(oecd_rwb$bic), "sd" = sd(oecd_rwb$bic), "min" = min(oecd_rwb$bic), "max" = max(oecd_rwb$bic), "source" = "CPDS"),
                                       cbind("var" = "Reg. unemployment rate", "mean" = mean(oecd_rwb$UNEM_RA), "sd" = sd(oecd_rwb$UNEM_RA), "min" = min(oecd_rwb$UNEM_RA), "max" = max(oecd_rwb$UNEM_RA), "source" = "OECD RWB"),
                                       cbind("var" = "Reg. share second. educ", "mean" = mean(oecd_rwb$EDU38_SH), "sd" = sd(oecd_rwb$EDU38_SH), "min" = min(oecd_rwb$EDU38_SH), "max" = max(oecd_rwb$EDU38_SH), "source" = "OECD RWB"),
                                       cbind("var" = "Reg. turnout", "mean" = mean(oecd_rwb$VOTERS_SH), "sd" = sd(oecd_rwb$VOTERS_SH), "min" = min(oecd_rwb$VOTERS_SH), "max" = max(oecd_rwb$VOTERS_SH), "source" = "OECD RWB")))

sum_stat_covars[ , 2:5 := lapply(.SD, function(x) round(as.numeric(as.character(x)),2) ), .SDcols= 2:5 ]

### SI Table C.2
print(xtable::xtable(sum_stat_covars, include.colnames=T, caption="Covariate Summary Statistics"), include.rownames=F)


                                    

### -------------
### Plot distribution 
### -------------

### Figure 1
pdf("out/fig1.pdf", width=5.5, height=3.5)  
(plot <- (ggplot(data=oecd_rwb, aes(y=diff_gini_idd, x=gini_before_tax, shape=elect_rules, fill=elect_rules, colour=elect_rules))
          + geom_jitter(alpha=0.5)
          + stat_smooth(method = "lm")
          + theme_bw() 
          + ylab("National-level redistribution")
          + xlab("Regional pre-tax gini coefficient")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), 
                  panel.grid = element_blank())
          + scale_shape_manual("Electoral regime", values=c(22,21))
          + scale_colour_brewer("Electoral regime", palette="Set1")
          + scale_fill_brewer("Electoral regime", palette="Set1")
))
dev.off()



### SI Figure A1
oecd_rwb_sub <- oecd_rwb[ ,.(country,gini_before_tax) ] 
oecd_rwb_sub$country <- factor(oecd_rwb_sub$country, levels=rev(unique(oecd_rwb_sub$country)))

pdf("out/si_fig_a1.pdf", width=5, height=5)  
(plot <- (ggplot(data=oecd_rwb_sub, aes(y=gini_before_tax, x=country))
          + geom_boxplot(fill="blue", outlier.shape = NA)
          + geom_point(size=1.5, shape=21, fill="gray", alpha=0.75)
          + theme_bw()
          + coord_flip()
          + ylab("Pre-tax and transfer gini")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), 
                  panel.grid = element_blank(), axis.title.y = element_blank())
))
dev.off()






### -------------
### Models 
### -------------

## OLS Models 
m1.1 <- lm(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu, oecd_rwb[  ])

m1.2 <- lm(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
             gov_left1 + coalition_gov, oecd_rwb[  ])

m1.3 <- lm(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
             gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, oecd_rwb[  ])

m1.4 <- lm(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
             gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ]  )

### IV models
m2.1 <- ivreg(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH
              | year_elec*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, data = oecd_rwb)

m2.2 <- ivreg(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH
              | year_elec*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, data = oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ]  )

m2.3 <- ivreg(diff_gini_idd ~ plurality_rule*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH
              | year_elec_grp*gini_before_tax + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, data = oecd_rwb  )

### First stage 
m2.1_FS <- lm(plurality_rule ~ year_elec + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov, data = oecd_rwb)

m2.2_FS <- lm(plurality_rule ~ year_elec + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, data = oecd_rwb)

m2.3_FS <- lm(plurality_rule ~ year_elec + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, data = oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ]  )

m2.4_FS <- lm(plurality_rule ~ year_elec_grp + gdp_log + unemp_rate + trade_union_density + eu + fed + pres_bin + bic + 
                gov_left1 + coalition_gov + UNEM_RA + EDU38_SH + VOTERS_SH, data = oecd_rwb  )


### Bootstrap cluster SE 
m1.1_formula <- paste(m1.1$call)[2]
m1.2_formula <- paste(m1.2$call)[2]
m1.3_formula <- paste(m1.3$call)[2]
m1.4_formula <- paste(m1.4$call)[2] 

m2.1_formula <- paste(m2.1$call)[2]
m2.2_formula <- paste(m2.2$call)[2]
m2.3_formula <- paste(m2.3$call)[2] 

m2.1_FS_formula <- paste(m2.1_FS$call)[2]
m2.2_FS_formula <- paste(m2.2_FS$call)[2] 
m2.3_FS_formula <- paste(m2.3_FS$call)[2] 
m2.4_FS_formula <- paste(m2.4_FS$call)[2] 

ols_mod_fun <- function(data,b,formula){
  d <- data[b,]  
  fit <- lm(formula, data=d)
  return(fit$coefficients)
}

iv_mod_fun <- function(data,b,formula){
  d <- data[b,] # 
  fit <- ivreg(formula, data=d)
  return(fit$coefficients)
}

set.seed(02138)
boot_ols_m1.1 <- boot(data=oecd_rwb, statistic=ols_mod_fun, R=2000, formula=m1.1_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")
boot_ols_m1.2 <- boot(data=oecd_rwb, statistic=ols_mod_fun, R=2000, formula=m1.2_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")
boot_ols_m1.3 <- boot(data=oecd_rwb, statistic=ols_mod_fun, R=2000, formula=m1.3_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")
boot_ols_m1.4 <- boot(data=oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ], statistic=ols_mod_fun, R=2000, formula=m1.4_formula, strata = as.factor(oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ]$ISO3), parallel = "multicore")

boot_iv_m2.1 <- boot(data=oecd_rwb, statistic=iv_mod_fun, R=2000, formula=m2.1_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")
boot_iv_m2.2 <- boot(data=oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ], statistic=iv_mod_fun, R=2000, formula=m2.2_formula, strata = as.factor(oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ]$ISO3), parallel = "multicore")
boot_iv_m2.3 <- boot(data=oecd_rwb, statistic=iv_mod_fun, R=2000, formula=m2.3_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")

boot_ols_m2.1_FS <- boot(data=oecd_rwb, statistic=ols_mod_fun, R=2000, formula=m2.1_FS_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")
boot_ols_m2.2_FS <- boot(data=oecd_rwb, statistic=ols_mod_fun, R=2000, formula=m2.2_FS_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")
boot_ols_m2.3_FS <- boot(data=oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ], statistic=ols_mod_fun, R=2000, formula=m2.3_FS_formula, strata = as.factor(oecd_rwb[ !ISO3 %in% c("SVN","NZL","IRL") ]$ISO3), parallel = "multicore")
boot_ols_m2.4_FS <- boot(data=oecd_rwb, statistic=ols_mod_fun, R=2000, formula=m2.4_FS_formula, strata = as.factor(oecd_rwb$ISO3), parallel = "multicore")

gc()


### Table 1  
stargazer(m1.1,m1.2,m1.3,
          m2.1,m1.4,m2.2,
          se = list(apply(boot_ols_m1.1$t,2,sd),
                    apply(boot_ols_m1.2$t,2,sd),
                    apply(boot_ols_m1.3$t,2,sd),
                    apply(boot_iv_m2.1$t,2,sd),
                    apply(boot_ols_m1.4$t,2,sd),
                    apply(boot_iv_m2.2$t,2,sd)),  
          digits.extra=0, digits=3, align=T, no.space=T, omit.stat=c("ser","F"), type="text",
          title = "Effect of Spatial Concentration of Inequality on Redistribution", table.placement = "h!",
          keep = c("plurality_rule", "gini_before_tax"),
         covariate.labels = c("Plurality rule","Pre-tax gini","Plurality rule $\\times$ Pre-tax gini"),
          dep.var.labels = "National redistribution",  
          add.lines = list(c("Mean DV", 
                             round(mean(m1.1$model$diff_gini_idd),3), 
                             round(mean(m1.2$model$diff_gini_idd),3), 
                             round(mean(m1.3$model$diff_gini_idd),3), 
                             round(mean(m2.1$model$diff_gini_idd),3),
                             round(mean(m1.4$model$diff_gini_idd),3), 
                             round(mean(m2.2$model$diff_gini_idd),3)),
                           c("No. countries", "\\multicolumn{1}{c}{25}","\\multicolumn{1}{c}{25}","\\multicolumn{1}{c}{25}",
                             "\\multicolumn{1}{c}{25}","\\multicolumn{1}{c}{22}","\\multicolumn{1}{c}{22}"),
                           c("F-Statistic","-","-","-",
                             round(summary(m2.1, diagnostics=T)$diagnostics[2,3],2),"-",
                             round(summary(m2.2, diagnostics=T)$diagnostics[2,3],2)),
                           c("Wu-Hausman","-","-","-",
                             round(summary(m2.1, diagnostics=T)$diagnostics[3,3],2),"-",
                             round(summary(m2.2, diagnostics=T)$diagnostics[3,3],2))))


### SI Table A3
covar.labFS <- c("Origin year electoral rule", 
                 "Log GDP","Unemployment rate","Union density","EU Member",
                 "Weak federalism","Strong federalism","Presidential system","Bicamerialism index",
                 "Share cabinet seats left parties","Coalition government",
                 "Regional unemployment rate","Regional share of workers with secondary education or higher","Regional turnout")

stargazer(m2.1_FS, m2.2_FS, m2.3_FS,
          se = list(apply(boot_ols_m2.1_FS$t,2,sd),
                    apply(boot_ols_m2.2_FS$t,2,sd),
                    apply(boot_ols_m2.3_FS$t,2,sd)),  
          digits.extra=0, digits=3, align=T, no.space=T, omit.stat=c("ser","F"), type="text",
          title = "IV First Stage Results: Effect of Year of Electoral Regime Origin on Electoral Regime Type", table.placement = "h!",
          covariate.labels = covar.labFS,
          dep.var.labels = "Plurality rule", 
          label = "tab:oecd_reg_mod1_IV_FS_v2",
          add.lines = list(c("No. countries", "\\multicolumn{1}{c}{25}", "\\multicolumn{1}{c}{25}", "\\multicolumn{1}{c}{22}")))


### SI Table A4
covar.labFS_v2 <- c("Origin electoral rule: 1900s", "Origin electoral rule: 1920s", "Origin electoral rule: 1940s", 
                    "Origin electoral rule: 1960s", "Origin electoral rule: 1980s", "Origin electoral rule: 2000s", 
                    "Log GDP","Unemployment rate","Union density","EU Member",
                    "Weak federalism","Strong federalism","Presidential system","Bicamerialism index",
                    "Share cabinet seats left parties","Coalition government",
                    "Regional unemployment rate","Regional share of workers with secondary education or higher","Regional turnout")

stargazer(m2.4_FS, se = list(apply(boot_ols_m2.4_FS$t,2,sd)),  
          digits.extra=0, digits=3, align=T, no.space=T, omit.stat=c("ser","F"), type="text",
          title = "IV First Stage Results: Effect of Year of Electoral Regime Origin on Electoral Regime Type", table.placement = "h!",
          covariate.labels = covar.labFS_v2,
          dep.var.labels = "Plurality rule",  
          add.lines = list(c("No. countries", "\\multicolumn{1}{c}{25}")))



### SI Table A5
covar.lab <- c("Plurality rule","Pre-tax gini", 
               "Log GDP","Unemployment rate","Union density","EU Member",
               "Weak federalism","Strong federalism","Presidential system","Bicamerialism index",
               "Share cabinet seats left parties","Coalition government",
               "Regional unemployment rate","Regional share of workers with secondary education or higher","Regional turnout",
               "Plurality rule $\\times$ Pre-tax gini")

stargazer(m1.1,m1.2,m1.3,
          m2.1,m1.4,m2.2, m2.3, 
          se = list(apply(boot_ols_m1.1$t,2,sd),
                    apply(boot_ols_m1.2$t,2,sd),
                    apply(boot_ols_m1.3$t,2,sd),
                    apply(boot_iv_m2.1$t,2,sd),
                    apply(boot_ols_m1.4$t,2,sd),
                    apply(boot_iv_m2.2$t,2,sd),
                    apply(boot_iv_m2.3$t,2,sd)), 
          digits.extra=0, digits=3, align=T, no.space=T, omit.stat=c("ser","F"), type="text",
          title = "Effect of Spatial Concentration of Inequality on Redistribution", table.placement = "h!",
          covariate.labels = covar.lab,
          dep.var.labels = "National redistribution",  
          add.lines = list(c("Mean DV", 
                             round(mean(m1.1$model$diff_gini_idd),3), 
                             round(mean(m1.2$model$diff_gini_idd),3), 
                             round(mean(m1.3$model$diff_gini_idd),3), 
                             round(mean(m2.1$model$diff_gini_idd),3),
                             round(mean(m1.4$model$diff_gini_idd),3), 
                             round(mean(m2.2$model$diff_gini_idd),3),
                             round(mean(m2.3$model$diff_gini_idd),3)),
                           c("No. countries", "\\multicolumn{1}{c}{25}","\\multicolumn{1}{c}{25}","\\multicolumn{1}{c}{25}",
                             "\\multicolumn{1}{c}{25}","\\multicolumn{1}{c}{22}","\\multicolumn{1}{c}{22}"),
                           c("F-Statistic","-","-","-",
                             round(summary(m2.1, diagnostics=T)$diagnostics[2,3],2),"-",
                             round(summary(m2.2, diagnostics=T)$diagnostics[2,3],2),
                             round(summary(m2.3, diagnostics=T)$diagnostics[2,3],2)),
                           c("Wu-Hausman","-","-","-",
                             round(summary(m2.1, diagnostics=T)$diagnostics[3,3],2),"-",
                             round(summary(m2.2, diagnostics=T)$diagnostics[3,3],2),
                             round(summary(m2.3, diagnostics=T)$diagnostics[3,3],2))))



### ----------------------
### Marginal effects plot
### ----------------------

var.level.m2.1 <- seq(min(oecd_rwb[ , gini_before_tax ], na.rm=T), 
                      max(oecd_rwb[ , gini_before_tax ], na.rm=T),
                      length.out = 25)

# Coefficients
estmean.m2.1 <- m2.1$coefficients
var.m2.1 <- vcov(boot_iv_m2.1) 

marg.eff_m2.1 <- m2.1$coefficients[2] + m2.1$coefficients[17]*var.level.m2.1

SEs_m2.1 <- rep(NA, length(var.level.m2.1))

for (i in 1:length(var.level.m2.1)){
  j <- var.level.m2.1[i]
  SEs_m2.1[i] <- msm::deltamethod (~  (x2) + (x17)*j , estmean.m2.1, var.m2.1, ses=T)
}

upper.m2.1 <- marg.eff_m2.1 + 1.96*(SEs_m2.1)
lower.m2.1 <- marg.eff_m2.1 - 1.96*(SEs_m2.1)

tab_m2.1 <- as.data.table(cbind(var.level=var.level.m2.1, lower = lower.m2.1, coef = marg.eff_m2.1, upper = upper.m2.1))

# Density plot
hist.out <- hist(oecd_rwb[ , gini_before_tax ], breaks=80, plot=F)
n.hist <- length(hist.out$mids)
dist <- hist.out$mids[2]-hist.out$mids[1]
hist.max <- max(hist.out$counts)
yrange <- c(tab_m2.1[ , upper], tab_m2.1[ ,lower])
maxdiff<-(max(yrange)-min(yrange))

# Create data frame
hist_m2.1 <- data.frame(ymin=rep(min(yrange)-maxdiff/5,n.hist),
                        ymax=hist.out$counts/hist.max*maxdiff/5+min(yrange)-maxdiff/5,
                        xmin=hist.out$mids-dist/2,
                        xmax=hist.out$mids+dist/2)

hist_m2.1 <- as.data.table(hist_m2.1)

pdf("out/fig2.pdf", width=4.5, height=3.5)  
(plot <- (ggplot()
          + geom_rect(data=hist_m2.1, aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax), fill="gray90", alpha=0.5, colour="gray70")
          + geom_hline(yintercept=0, colour = "gray50", linetype="dashed", size=.5)
          + geom_ribbon(data=tab_m2.1, aes(y=coef, x=var.level, ymin=lower, ymax=upper), alpha=0.3, color=NA, fill="blue")
          + geom_line(data=tab_m2.1, aes(y=coef, x=var.level), linewidth=0.9)
          + ylab("Marginal effect of plurality\nelectoral rule on redistribution")
          + xlab("Regional pre-tax inequality")
          + theme_bw()
          + theme(legend.position="bottom", panel.grid=element_blank(), panel.border = element_rect(color="gray"))
))
dev.off()


